#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <vector>



#define  HUBBLE      1.0     /* Hubble constant in 100km/sec/Mpc */  
#define  PI  3.1415926

#define  GRAVITY			6.672e-8
#define  SOLAR_MASS		1.989e33
#define  SOLAR_LUM		3.826e33
#define  RAD_CONST		7.565e-15
#define  AVOGADRO			6.0222e23
#define  BOLTZMANN		1.3806e-16
#define  GAS_CONST		8.31425e7
#define  C					2.9979e10
#define  PLANCK			6.6262e-27
#define  CM_PER_MPC		3.085678e24
#define  PROTONMASS		1.6726e-24
#define  ELECTRONMASS	9.10953e-28
#define  THOMPSON			6.65245e-25
#define  ELECTRONCHARGE	4.8032e-10

#define  SEC_PER_MEGAYEAR   3.155e13
#define  SEC_PER_YEAR       3.155e7


using namespace std;

double G;			// gravitational constant
double H0;			// Hubble constant
double UnitTime_in_s;
double UnitMass_in_g;
double UnitLength_in_cm;
double UnitVelocity_in_cm_per_s;
double UnitTime_in_Megayears;

double UnitPressure_in_cgs, UnitDensity_in_cgs, UnitCoolingRate_in_cgs, UnitEnergy_in_cgs;

double CC = 21.771771;      // halo concentration 
double V200 = 166.907314; // circular velocity v_200 (in km/sec)

double MD = 0.027746;
double MB = 0.0092511;
double LAMBDA=0.033;    // spin parameter  
double JD = 0.027746;
double DiskHeight = 0.2;
double BulgeSize = 0.2;

double R200;
double M200;
double RS;
double RH;	// scale length of Hernquist profile

double M_DISK;
double M_BULGE;
double M_HALO;

double HubbleParam;

double Z0; // sets disk thickness for stars 
double  A;  // sets bulge size 

void structure(void) {
  //double jhalo, jdisk, jd;
  //double hnew, dh;

  R200 = V200 / (10 * H0);
  M200 = pow(V200, 3) / (10 * G * H0);

  RS = R200 / CC;
  M_DISK = MD * M200;
  M_BULGE = MB * M200;

  //M_GAS = M_DISK * GasFraction;
  M_HALO = M200 - M_DISK - M_BULGE;

  RH = RS * sqrt(2 * (log(1 + CC) - CC / (1 + CC)));	// scale length of Hernquist profile 

  //double jhalo = LAMBDA * sqrt(G) * pow(M200, 1.5) * sqrt(2 * R200 / fc(CC));
  //double jdisk = JD * jhalo;

  //printf("fc=%g\n", fc(CC));

  //double halo_spinfactor = 1.5 * LAMBDA * sqrt(2 * CC / fc(CC)) * pow(log(1 + CC) - CC / (1 + CC), 1.5) / gc(CC);

  //printf("halo_spinfactor %g\n", halo_spinfactor);
  printf("R200: %g\n", R200);
  printf("M200: %g\n", M200);
  printf("RH: %g\n", RH);
  printf("Total mass is: %g\n", M200);
	/*

  H = sqrt(2.0) / 2.0 * LAMBDA / fc(CC) * R200;	// first guess for disk scale length 
*/
	
	H = 1.9;
  Z0 = DiskHeight * H;		// sets disk thickness for stars 
  A = BulgeSize * H;		// sets bulge size 
  
	/*
  if(M_DISK > 0)
    {
      do
	{
	  jd = disk_angmomentum();	// computes disk momentum
	  
	  hnew = jdisk / jd * H;
	  
	  dh = hnew - H;
	  
	  if(fabs(dh) > 0.5 * H)
	    {
	      dh = 0.5 * H * dh / fabs(dh);
	    }
	  else
	    dh = dh * 0.1;
	  
	  H = H + dh;
	  
	  printf("Jd/J=%g   hnew: %g  \n", jd / jhalo, H);
	  
	  
	  Z0 = DiskHeight * H;	// sets disk thickness 
	  A = BulgeSize * H;	// sets bulge size 
	}
      while(fabs(dh) / H > 1e-5);
    }
    */
}


void init_units(void) {

  UnitMass_in_g = 1e10 * SOLAR_MASS;	/* 10^10 solar masses */
  UnitLength_in_cm = CM_PER_MPC / 1000;	/* 1 kpc */
  UnitVelocity_in_cm_per_s = 1e5;	/* 1 km/sec */

  HubbleParam = 0.7;

  UnitTime_in_s = UnitLength_in_cm / UnitVelocity_in_cm_per_s;
  H0 = HUBBLE * 100 * 1e5 / CM_PER_MPC / UnitVelocity_in_cm_per_s * UnitLength_in_cm;
  G = GRAVITY / pow(UnitLength_in_cm, 3) * UnitMass_in_g * pow(UnitTime_in_s, 2);

  UnitTime_in_Megayears = UnitTime_in_s / SEC_PER_MEGAYEAR;

  UnitDensity_in_cgs = UnitMass_in_g / pow(UnitLength_in_cm, 3);
  UnitPressure_in_cgs = UnitMass_in_g / UnitLength_in_cm / pow(UnitTime_in_s, 2);
  UnitCoolingRate_in_cgs = UnitPressure_in_cgs / UnitTime_in_s;
  UnitEnergy_in_cgs = UnitMass_in_g * pow(UnitLength_in_cm, 2) / pow(UnitTime_in_s, 2);


  printf("%g %g %g \n", G, H0, UnitTime_in_Megayears);
}


double comp_Dphi_R_bulge(double RR, double zz)
{
  double m;
  double r;

  r = sqrt(RR * RR + zz * zz);

  m = M_BULGE * r * r / pow(A + r, 2);

  if(r > 0)
    return G * RR / (r * r * r) * m;
  else
    return 0;
}
/*

double halo_mass(double r)
{
  return M_HALO * pow(r / (r + RH), 2);
}

double comp_Dphi_R_halo(double R, double z)
{
  double r, M_r;
  double halo_mass(double r);


  r = sqrt(R * R + z * z);

  M_r = halo_mass(r);

  if(r > 0)
    return G * R / (r * r * r) * M_r;
  else
    return 0;
}


double comp_Dphi_R_disk_tree(double RR, double zz)
{
  double pos[3], acc[3];
  int nint;

  if(M_DISK > 0)
    {
      pos[0] = RR;
      pos[1] = 0;
      pos[2] = zz;
      
      nint = force_treeevaluate(pos, 0.01 * H, acc);
  
      return -G * acc[0];
    }
  else
    return 0;
}


double comp_Dphi_R(double R, double z) {
  return comp_Dphi_R_disk_tree(R, z) + comp_Dphi_R_halo(R, z) + comp_Dphi_R_bulge(R, z);
}

*/

double epicyclic_kappa2(double R) {
  double dR, dphi, dphi_;

  if(R > 0)
    {
      dR = R * 0.01;

      dphi = comp_Dphi_R(R, 0);

      dphi_ = comp_Dphi_R(R + dR, 0);

      return 3 * dphi / R + (dphi_ - dphi) / dR;

    }
  else
    return 0;
}



void plot_toomre_stability(FILE * fd)
{
  double *Q;
  int i;
  double Sigma0;
  int count;

  for(i = 0, count = 0; i <= RSIZE; i++)
    if(list_R[i] <= 6 * H)
      count++;

  Sigma0 = (M_DISK) / (2 * PI * H * H);
  Q = dvector(0, RSIZE);


  for(i = 0; i < count; i++)
    Q[i] = RadialDispersionFactor * sqrt(VelDispRz_disk[i][0]) * sqrt(epi_kappa2[i]) / (3.36 * G * Sigma0 * exp(-list_R[i] / H));


  fprintf(fd, "\n%d\n", count);

  for(i = 0; i < count; i++)
    fprintf(fd, "%g\n", list_R[i]);

  for(i = 0; i < count; i++)
    fprintf(fd, "%g\n", Q[i]);
}

*/
